Curvature-based interface restoration algorithm using phase-field equations

In this study, we propose a restoration algorithm for distorted objects using a curvature-driven flow. First, we capture the convex-hull shaped contour of the distorted object using the mean curvature flow. With the supplemented mass on the captured feature, we evolve the constraint mean curvature flow to a steady state, preserving the non-distorted region. With respect to the mass, we select a restorative shape by considering the square of the curvature derivative. The Allen–Cahn and Cahn–Hilliard equations are applied to the generated restored image from the implicit curvature motions represented by the order parameter. We impose the Dirichlet boundary condition for the order parameter and the Neumann boundary for the chemical potential to fix the feature and to inherit the mass conservation, respectively. We provided examples of the restoration of half-circle and parentheses-shaped objects to reconstruct a circle shape.


Introduction
Image data in electronic formats can be damaged in storage or transmission [1].In medical tomography, imaging scans can be easily distorted by artificial effects, causing interference and interruption [2].For example, in surgical imaging, a part of a medical image showing the skeleton could be accidentally damaged, which might cause difficulty in distinguishing the characteristic signs of some conditions or diseases [3].Therefore, considerable effort has been devoted to the study of the reconstruction of distorted images.
Image inpainting techniques reconstruct image data using information around damaged parts.shows an image reconstructed using an image inpainting technique.Many models have been proposed to obtain a reasonable reconstruction [4][5][6][7][8] based on modifications of the Allen-Cahn and Cahn-Hilliard equations [9][10][11].These successive binary image inpainting methods have been extended to a wide variety of applications [12][13][14][15] such as colored images and slice data in three-dimensional problems.
By contrast, an image domain (surface mesh) could be destroyed along with the image value, as shown in Fig 2. This has been investigated in terms of hole-filling or surface inpainting.In 3D cases, a volumetric reconstruction may be considered.However, in the present work, we discuss only surface restoration [16][17][18][19][20][21][22][23].We refer interested readers to a survey of the relevant literature [24].In [16], volumetric diffusion method was proposed as an algorithm to restore a complex surface.Moreover, several methods have been proposed that considered geometric properties such as curvature [17,18,21,23].In [19], a mesh generation and self-intersection from the projected boundary points on a plane were used to update surface mesh points inside a hole in a piecewise manner.In [20], an advancing front mesh algorithm was applied to fill in a hole, and the Poisson equation was used to heal over inside points.In [22], a sparsity-based hole-filling algorithm was proposed to find a low-dimensional intrinsic structure from a Laplacian eigenbasis dictionary.
In comparison, to distinguish surface restoration with image inpainting, the latter indicates the reconstructing of function values at grid points in the damaged region (in the case that we already know the locations of the grid points on a flatted 2D image), whereas the former reconstructs grid points for a damaged curved 2D image.To discuss this clearly, we consider the concept of interface restoration, which aims to recover a parametric curve from a destroyed curve.For example, given the open curve with two endpoints in Fig 3(a), we aim to find the closed curve that links the two points, where we assume that the original is closed, smooth, and non-twisted.
In this study, we propose a curvature-based interface restoration method.The desired curve is selected by the minimized curve of the curvature variation.We provide the proposed algorithm with a simple example in which a circle shape is partially lost.For example, refer to  had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
V 0 .Moreover, we define D ΔV as the curvature minimization shape with the area of V 0 + ΔV.Then, we select the desired curve by considering the curvature variation with respect to the increment of volume ΔV.
We employ the Allen-Cahn (AC) equation ϕ t = −μ and the Cahn-Hilliard (CH) equation ϕ t = Δμ, where μ = ϕ 3 − ϕ − � 2 Δϕ is the chemical potential.These equations are considered to use the properties of curvature-driven flows [25][26][27], which were originally introduced to describe motion of interfaces in materials [28,29].Various studies are being conducted on the AC equation that conserves mass [30].This can be used instead of the CH equation, but since the characteristics of the two equations are slightly different [31], so the CH equation is used in this research.In addition, we define a narrow band shape to represent the given damaged curve and the outer boundary to approximate the desired restored curve.We consider the Dirichlet boundary condition for both the AC and CH equations to fix the boundary near the undamaged curve, and the Neumann boundary condition for the CH equation to preserve the supplemented area.We described the detailed algorithm in the Numerical method section.
The remainder of this study is organized as follows.In the Numerical method section, we present the proposed numerical method with restoration algorithms defining capture and restorations step.In the Numerical results section, we describe the results of numerical simulations including half-circle and parentheses shapes.Finally, Conclusion section presents a discussion of our conclusions.

Numerical method
In this section, we describe each step in interface restoration in detail.The overall process is shown in The proposed restoration algorithm consists of the following two sub-steps.The first step (capture) identifies the overall structure, including the missing part, and the second step (reconstruction) restores the shape through curvature minimization.
1.The capture step sets the order parameter �ðx; tÞ to capture a given damaged feature by using the AC equation with one Dirichlet boundary condition.

Capture step
The governing equation of the capture step is the AC equation with one Dirichlet boundary condition for �ðx; tÞ, as follows.
@�ðx; tÞ @t ¼ À ð�ðx; tÞÞ 3 À �ðx; tÞ where � is the coefficient related with an interface thickness.The initial condition �ðx; 0Þ is satisfied in that its values are defined as 1 for x 2 O d and its zero-contour covers O d .At this point, ( 3) and ( 4) can be discretized on G h as follows.
where � n ij is the approximations of �ðx ij ; nDtÞ, and Δt is the time step size.The capture step is complete when � n ij reaches a steady state.

Reconstruction step
The reconstruction step is governed by the types of the CH equations using the boundary control function Here, the initial condition is set by the steady-state solution of the capture step.
The discretization of the system ( 7) and ( 9) on G h can be expressed as follows. where Because the summation by parts is satisfied for m n ij [36], discrete mass conservation holds as follows.

Numerical results
We set the tolerance for the numerical simulation as tol = 0.5 � 10 −6 .

Grid-independent test
To confirm our method, we first perform the numerical simulation that gives independency of the number of grid points before full restoration steps.Here, we set a zero-Neumann boundary condition and the initial condition where dðx; yÞ ¼ 0:35 À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À 0:5Þ 2 þ ðy À 0:5Þ The slope of graph represents the second-order accuracy consistenting with the theory in [32].Moreover, it indicates a grid-independency that the scaling of error is smaller than 10 −1 in spite of the small number of the grid point.

Half-circle shape
We begins by evolving the solution of the AC equation with a zero-Neumann boundary condition and the initial condition same as above.Fig 7 shows the evolution of the solution of the AC equation in blue and the obstacle shape in red.Because the AC equation can be explained as the mean curvature flow, the initial circle rapidly shrinks before the iso-contour contacts the given obstacle.Then, they evolve into half-circle shape, preserving the convex hull property.
Next, using the steady-state solution c from the AC equation, we evolve the CH equation with the same obstacle to the steady state.For the numerical simulation of a desired volume V β with a parameter of β, we set an initial condition where dðx; yÞ ¼ ðr À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À 0:5Þ 2 þ ðy À 0:5 À r À 4DyÞ 2 q , r ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi bðV b þ V 0 Þ=ð2pÞ q , and To investigate the curvature distribution with respect to the volume, we compute the curvature on the zero-contour of the steady-state solution with various β.To numerically compute the curvature, we generate the contour points using the DistMesh algorithm [37].Fig 9 shows the curvature plots for the various values of β.The blue line shows the zero-contour of the final solution at t = t f , and the red line represents the curvature value on the zero-contour.With increasing β, the distribution of curvature is smeared and the high magnitude near the endpoints decrease.
To characterize the steady-state shape using the curvature distribution, we consider the integration of the square of curvature derivative as  condition.Therefore, a standard is required to make a reasonable choice, and we thus adopt the steady-state solution for β = 0.9 and β = 1.

Parentheses shape
We begin by evolving the solution of the AC equation with a zero-Neumann boundary condition and the initial conditions, as given below.where dðx; yÞ ¼ 0:35 À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À 0:5Þ The obstacle is defined as the middle of parentheses shapes describing parts of a circle with radii of 0.2 and 0.24.The center of the circles is the center of the computational domain.
Fig 11 shows the evolution of the solution for the AC equation in blue and the obstacle shape in red.Similarly, for the results of the previous section, the initial circle rapidly shrinks before the iso-contour contacts the given obstacle and then evolves into the tablet shape, preserving the convex hull property.
Next, using the steady-state solution c from the AC equation, we evolve the CH equation with the same obstacle to the steady state.For the numerical simulation of a desired volume V β with a parameter of β, we set an initial condition where ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À 0:5Þ 2 þ ðy À 0:7 À 0:5r À 4DyÞ 2 q ; ð19Þ d 2 ðx; yÞ ¼ r À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À 0:5Þ 2 þ ðy À 0:3 þ 0:5r þ 4DyÞ 2 q ; ð20Þ r ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi bðV b þ V 0 Þ=ð4pÞ q , and V 0 ¼ meanð cÞ.With the initial condition (18), we consider the steady-state solution of the constraint mean-curvature flow with the specific volume V β .

Skull shape in CT image
The previous research [38], while conducting research on the development of skull protheses through an active contour model, restored the part of the skull if it takes damaged (missing) part from the corresponding CT image (skull).To confirm our proposed algorithm, we also choose the same CT image (Fig 15

Conclusion
In this study, we proposed a restoration algorithm for a distorted object using curvaturedriven flows.We captured the outer contour using the Allen-Cahn equation with the one-Dirichlet boundary condition for the distorted object.Subsequently, we the Cahn- Hilliard equation to a steady state with the supplemented mass on the captured feature.Imposing the one-Dirichlet boundary condition for the order parameter and zero-Neumann boundary condition for the chemical potential, we fixed the feature and inherited the mass conservation.Finally, we chose the restored shape by considering the square of the curvature derivative.We provided two restoration examples for half-circle and parentheses objects showing how the proposed approach reconstructed a circle shape.The numerical results showed that the curvature distribution became more diffused with increasing supplemented mass, and it converged to the circle shape.The approximated retoration gave expected results; however, it may be broken when we set the too large mass since the one-Dirichlet boundary condition cannot exactly capture the fixed endpoint.We leave the problem of fixing the endpoints as a subject for future work.
Fig 1 shows a well-known 2D schematic of a conventional method to generate gray values in part of a vacant region.Fig 1(a) shows a given image containing a missing domain, and Fig 1(b)

Fig 3 .
First, we denote the domain generated by connecting two points as D 0 , and its area by

Fig 4 .
First, we consider a narrow band to approximate the given shape in Fig 4(a).The steady-state solution of the AC equation is the shape of D 0 in Fig 4(b), and we calculate its area V 0 , which is referred to as the capture step.Evolving the CH equation with the supplemented area V 0 + ΔV, we obtain a steady-state solution D ΔV , referred to as the reconstruction step in Fig 4(c).Finally, we select a reasonable shape D DV C depending on the curvature derivation in Fig 4(d).

Fig 3 .Fig 4 .Fig 5 .
Fig 3. Interface restoration; (a) open curve, (b) closed region D 0 obtained by filling the boundary with an area of V 0 , (c) curvature minimization shape D ΔV with an area of V 0 + ΔV, (d) selected result with a reasonable area of V 0 + ΔV C .https://doi.org/10.1371/journal.pone.0295527.g003 Moreover, we extend the definition of G ij to the cell-edged grid points as follows (See Fig 5 (c)).

Fig 6 .Fig 7 .
Fig 6.Relative l 2 -error with respect to spatial step size.https://doi.org/10.1371/journal.pone.0295527.g006 The corresponding values of D k of Fig 9 are 10900, 5190, 1750, 330, 8.85, and 1.14, respectively.Fig 10 shows the curvature variation D k with the various β, the value of which monotonically decreases.It should be noted that choice of a too large area V 0 (1 + β) could make an interface diverge from the target shape in spice of the obstacle and the Dirichlet boundary

Fig 9 .Fig 10 .
Fig 9. Curvature values on the zero-contour at the final time.The corresponding values of D k : (a) 10900, (b) 5190, (c) 1750, (d) 330, (e) 8.85, and (f) 1.14.https://doi.org/10.1371/journal.pone.0295527.g009 Fig 12 shows the evolution of the solution for the CH equation with several parameters for β.The last column indicates the numerical steady-state solution such that the consecutive error of the solution is less than tol.The red-solid and black-dashed curves represent the obstacle and the final solution c for the AC equation, respectively.Owing to the nature of the mean curvature flow, the iso-contour of the solution for the CH equation near the obstacle tends to remain, and the other regions are searched for the optimized circle-like shape.By increasing the value of β, the final solutions converge to the circle.Fig 13 shows the curvature plots for the selected parameters of β.The blue line is the zerocontour of the final solution at t = t f and the red line presents the curvature value.The corresponding values of D k are 6350, 752, 14.5, and 1.30, respectively.As shown in Fig 14, the curvature variation D k monotonically decreases with the various β, as shown in as Fig 10.In this case, the steady-state solutions for β = 0.35 and β = 0.4 appear to be the reasonable choices.

Fig 11 .
Fig 11.Evolution of the zero-contour for AC equation.https://doi.org/10.1371/journal.pone.0295527.g011 (a)) and perform the restoration with two different damaged scenarios: front left cut out [case 1] in Fig 15(b) and back of the head cut out [case 2] in Fig 15 (c).The numerical results are shown in Fig 15(d)-15(f).Note that both results have good agreements with the original skull shape.We upload the basic data and visualization code for this simulation to the corresponding author's web page (https://sites.google.com/view/yhchoi/code)and DRYAD site (DOI: 10.5061/dryad.fbg79cp2g).